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S Abstract. We study effects of turbulent mixing on the critical behaviour of a 
nonequilibrium system near its second-order phase transition between the absorbing 
and fluctuating states. The model describes the spreading of an agent (e.g., infectious 
disease) in a reaction-diffusion system and belongs to the universality class of the 
directed bond percolation process, also known as simple epidemic process, and is 
equivalent to the Reggeon field theory. The turbulent advecting velocity field is 
modelled by the Obukhov-Kraichnan's rapid-change ensemble: Gaussian statistics 
with the correlation function (vv) oc 6(t — t') k~ d ~^ , where k is the wave number 
£H . and < £ < 2 is a free parameter. Using the field theoretic renormalization group 

we show that, depending on the relation between the exponent £ and the spatial 
dimension d, the system reveals different types of large-scale asymptotic behaviour, 
associated with four possible fixed points of the renormalization group equations. In 
addition to known regimes (ordinary diffusion, ordinary directed percolation process, 
and passively advected scalar field), existence of a new nonequilibrium universality 
class is established, and the corresponding critical dimensions are calculated to first 
' order of the double expansion in £ and e = 4 — d (one-loop approximation). It turns 

out, however, that the most realistic values £ = 4/3 (Kolmogorov's fully developed 
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turbulence) and d — 2 or 3 correspond to the case of passive scalar field, when 
■ the nonlinearity of the Reggeon model is irrelevant and the spreading of the agent 

, is completely determined by the turbulent transfer. 
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1. Introduction 

Over the past few decades, constant interest has been attracted by the spreading 
processes and corresponding nonequilibrium phase transitions; see Refs. [T] — [10] and 
the literature cited therein. Spreading processes are ubiquitous in the Nature and 
are encountered in physical, chemical, biological, ecological and sociological systems: 
autocatalytic reactions, percolation in porous media, forest fires, epidemic diseases, and 
so on. For definiteness, in the following we will use the terminology of the latter case 
(spreading of a disease). 

Depending on the conditions, the spreading of an agent (disease) can either continue 
over the whole population, or terminate after some time. In the first case, the system 
evolves to a stationary (but not thermally equilibrium) active state, in which the sick 
and healthy individuals coexist, passing repeatedly through infection and healing events, 
and their densities are fluctuating (random) quantities. In the second case, when 
the probability of being healed is large enough in comparison with the probability 
to be infected, the system is "trapped" in an absorbing (inactive) state where all 
the individuals are healthy and all the fluctuations cease completely. The transitions 
between these fluctuating and absorbing phases are continuous; they are especially 
interesting as examples of nonequilibrium critical behaviour. 

By analogy with the equilibrium second-order transitions pJ],[T2] it is expected that, 
near the critical point, many details of a specific spreading process become irrelevant, so 
that the critical behaviour of different systems appears identical and can be described 
by a certain universality class. The aim of the theory is to identify possible universality 
classes and to calculate their universal characteristics (such as critical exponents, scaling 
functions) on the basis of an appropriate theoretic model and within a controlled 
approximation or a regular perturbation scheme. 

As a rule, the spreading phenomena are modelled by various stochastic reaction- 
diffusion processes on a lattice; see Ref. [1] for a detailed discussion. In the continuum 
limit, they can be mapped onto certain field theoretic models with the aid of special 
techniques [131 EE] • Then the powerful tools [HI [12] of the field theoretic renormalization 
group (RG) can be applied to the investigation of their critical behaviour. 

The most typical processes belong to the universality class of the so-called directed 
bond percolation process PQ [3], which in the field theoretic formulation is equivalent 
to the well-known Reggeon field theory [15]. It was conjectured in [J| that the critical 
behaviour always belongs to this "DP class" provided the absorbing state of the process 
is unique, the order parameter is one-dimensional, there are no specific symmetries, 
no coupling with additional "slow" degrees of freedom and no long-range interactions. 
Thus the DP process is expected jl] to play in the theory of noneqilibrium phase 
transitions the same paradigmatic role as the standard A0 4 model [HI [12] does in the 
theory of equilibrium critical behaviour. It is also sometimes referred to as simple 
epidemic process with recovery or as Gribov's process; the stochastic version of the 
Schlogl's first reaction also belongs to this class [T] — [3] . The corresponding critical 
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behaviour is rather well understood: above the upper critical dimension d c = 4 the 
critical exponents are given by the mean-field theory, and below d c they have been 
calculated to second order of the expansion in the deviation e = d c — d of the spatial 
dimension d from the upper critical value; see the discussion and references in Ref. [3]. 
Recently, the first experimental observation of this universality class was achieved in the 
transition between two topologically different turbulent states of electrohydrodynamic 
convection of a nematic liquid crystal [5], with a firm agreement between theoretical 
and experimental values of all the critical exponents. 

Extensive numerical and analytical investigations of various models of the spreading 
processes have shown that effects of immunization pQ [3], [7], long-range interactions [6] 
[8], presence of quenched (static) disorder j9] or critical fluctuations [10] in the medium 
where the spreading occurs can lead to appearance of completely new universality 
classes, with different values of the upper critical dimensions d c , different sets of critical 
exponents and interesting crossovers between various critical regimes. 

In this paper we will consider the spreading of a nonconserved agent in a turbulent 
medium and study the effects of turbulent stirring and mixing on the critical behaviour 
near the phase transition between the absorbing and the fluctuating phases. Turbulence 
is hardly avoidable in the chemical catalytic reactions and can play decisive role if a 
disease is transferred by flying insects or birds. 

Much work has been devoted to the investigation of the effects of various kinds 
of imposed deterministic or chaotic flows (laminar shear flows, turbulent convection 
and so on) on the behaviour of the critical or nearly critical fluids, like binary liquid 
mixtures near and below their consolution points; see e.g. the papers PS] [25] and 
references therein. It was shown that the flow can destroy the usual critical behaviour, 
typical of the A0 4 model, which changes to the mean-field behaviour [T61 [T7] or to a 
complex behaviour described by new nonequilibrium universality classes [T9l [20| [21]. 
This problem is closely related to another interesting issue: the effects of imposed flows 
on the dynamics of phase ordering - the growth of order through domain coarsening 
(spinodal decomposition) and formation of various nontrivial spatiotemporal patterns; 
see [22] [25]. 

Theoretical description of fully developed turbulence on the basis of a microscopic 
dynamic model (like e.g. the stirred Navier-Stokes equation) remains essentially an 
open problem [26]. In this paper, we will describe the turbulent mixing by a Gaussian 
velocity ensemble with prescribed statistics, which dates back to Obukhov, Batchelor, 
Kazantsev and Kraichnan and now is widely known as the rapid-change model |27] . In 
spite of their relative simplicity, the models of passive (no feedback on the velocity) 
scalar fields advected by such "synthetic" velocity ensembles have attracted serious 
attention recently because of the deep insight they offer into the origin of intermittency 
and anomalous ( mult i) scaling of the genuine turbulent heat or mass transport; see the 
review paper [27] and references therein. The RG approach to that problem is reviewed 
in Ref. [2E]. 

In the rapid-change model, the pair velocity correlation function is taken in the 
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form (vv) oc S(t — t') k~ d ~^, where k is the wave number and < £ < 2 a free parameter 
with the most realistic (Kolmogorov) value £ = 4/3. Vanishing of the correlation time 
ensures the Galilean symmetry of the problem, while a power-law dependence on the 
wave number k mimics the real scaling properties of fully developed turbulence. 

The plan of the paper is as follows. In section [2] we present the detailed description 
of the model and its field theoretic formulation. In section [3] we analyze the ultraviolet 
(UV) divergences of the model, relaying upon the canonical dimensions and additional 
symmetry considerations. We show that the model is multiplicatively renormalizable 
and write the renormalized action functional. 

In sectionHJwe derive the RG equations and introduce the RG functions {(3 functions 
and anomalous dimensions 7). In section [5] we identify four possible infrared (IR) 
attractive fixed points of the RG equations, and identify their ranges of stability in the 
£-£ plane. 

These fixed points correspond to four possible critical regimes of the model, with 
different sets of critical dimensions, as discussed in section [6] Three of them correspond 
to already known regimes: free (Gaussian) field theory, linear passive scalar advection 
by the Obukhov-Kraichnan ensemble (the nonlinearity in the agent's density appears 
irrelevant), and the ordinary Gribov process (mixing by the turbulent field is irrelevant). 
The fourth fixed point corresponds to a new nonequilibrium universality class, with the 
new set of critical exponents that depend on the both parameters e and £ and can be 
systematically calculated as double series in these parameters. 

The practical calculation of the renormalization constants, RG functions, regions 
of stability and critical dimensions is accomplished to the leading-order (one-loop) 
approximation; some of the results, however, are exact (valid to all orders of the double 
e-£ expansion). 

In section [6] we also consider, as a special consequence of the general scaling 
relations, the temporal evolution of a cloud of the advected agent's particles (infected 
individuals), which differs, in general, from the well-known "1/2 law" for ordinary 
diffusion. Section [7] is reserved for discussion and conclusion. 

2. Description of the model. Field theoretic formulation 

In the Lagrangian formulation, spreading of an agent is described by a continuous 
stochastic diffusion-reaction equation of the form 

d t iP(t, x) = A { (-r + <9 2 Mt, x) - g ^ 2 (t, x)/2} + ((f x), (2.1) 

where ip(t, x) > is the agent's density, d 2 is the Laplace operator, Ao and go are positive 
parameters, tq — deviation of the infection probability from its critical value (analog 
of To oc (T — T c ) in equilibrium systems). The random Gaussian noise C(A X ) with zero 
mean and a given correlation function 

(C(t, x)C(t', x')> = g X ^(f x)5(t - t')5^( X - x') (2.2) 

mimics fluctuations in the system; al is the dimension of the x space. 
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According to the general theorem [T3] (see also the monographs [HI [12]), stochastic 
problem (12.11) . (12.21) is equivalent to the field theoretic model of the doubled set of fields 
with action functional 

Sty, tf) = V f (-a + V? 2 - A r )^ + ^ ((^) 2 ^ - V'V 2 ) • (2.3) 

Here, tp^ = ip'(t,x) is the auxiliary "response field" and the integration over the 
arguments of the fields is implied, for example 

^d t tp = J dt J d-x^(t,x)d t ip(t,x). 

This means that statistical averages of random quantities in the original stochastic 
problem 02. ip . (I2.2I) can be represented as functional averages over the full set of fields 
with the weight exp5($), and can therefore be viewed as the Green functions of the 
field theoretic model with the action ( 12.3I) . In particular, the linear response function of 
the problem (I2.ll) . ( 12. 2ft is given by the Green function G = (ip^ip) of the model ( 12.31) . 

More rigorous derivation of the field theory is based on the so-called master 
equation for the original reaction-diffusion process on the lattice, which involves a 
microscopic Hamiltonian written in terms of the creation-annihilation operators with 
further representation of the corresponding evolution operator by the coherent-state 
functional integral [14J. For the case at hand the resulting action functional coincides 
with ( 12. 31) up to irrelevant terms; see also discussion of that approach in Refs. [2], [3]. 

One can argue that in the perturbation theory the condition ip > can be neglected 
[1]; then the model (I2.3P becomes equivalent to the Reggeon field theory [15] and acquires 
the symmetry with respect to the transformation 

V(t,x) -+V f (-t,-x), ^(t,x) ->^(-t,- x ), go^-go- (2.4) 

Reflection of the constant go is i n fact unimportant because, as can easily be seen, the 
actual expansion parameter in the perturbation theory is g% rather than g Q itself. 

The model (12.31) corresponds to a standard Feynman diagrammatic technique with 
the only bare propagator Go = (ipip')o and the two triple vertices ~ (ip^) 2 ^, tp^ifj 2 In 
the time-momentum and frequency-momentum representation Go has the forms 

G Q (t, k) = 9(t) exp {-A (A; 2 + r )} <- G (u, k) = , 1 2 (2.5) 

— ILO + Aq\K + To) 

Here #(...) is the Heaviside step function, so that the propagator (12. 5p is retarded. 
Then from the analysis of the diagrams one can check that the Green functions built 
solely from the field ip or solely from if)^ necessarily contain closed circuits of retarded 
propagators (12.51) and therefore vanish identically. For the functions (if)' . . . ip ') this fact 
is a general consequence of the causality, which is valid for any stochastic model; see 
e.g. the discussion in [12] . Then vanishing of the functions (tp . . .ip) can be viewed as a 
consequence of the symmetry (12.41) . 

The stability in dynamical models like ( 12.31) implies that all the small fluctuations 
are damped out, so that the exact response function G = (if)^i/>) must decay for t — > oo; 
see e.g. sec. 5.5 in [12]. Then from the expression ( 12.51) . which is the zero order 
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approximation for G, we conclude that the stability of the perturbative stationary state 
is lost for Tq = 0, and for r < the perturbations with smallest momenta grow in time. 
This growth is stabilized by the appearance of the nonzero constant mean (if)) and the 
higher-order correlation functions of the agent field if>, so that the symmetry (12.41) is 
spontaneously broken. This is exactly the phase transition from the absorbing (normal) 
to the fluctuating (anomalous) states. 

Coupling with the velocity field v = {i>j(i, x)} is introduced by the replacement 

d t ^Vt = d t + v i d i , (2.6) 

in (12.11) and (12. 3ft . where d{ = d/dxi and Vt is the Lagrangian (Galilean covariant) 
derivative. We will consider the case of incompressible flow, then the velocity field is 
divergence-free (transverse): diVi = 0. In the real problem, the field v(i, x) satisfies the 
Navier-Stokes equation. We will employ the rapid-change model, where the velocity 
obeys a Gaussian distribution with zero mean and correlation function 

(Vi(t, x)Vj(t', x')) = 5(t - t') Dij(T), r = x - x' (2.7) 
alk 

where i^j(k) = 6{j — kikj/k 2 is the transverse projector, D > is an amplitude factor 
and < £ < 2 is a free parameter with the most realistic ( "Kolmogorov" ) value £ = 4/3. 
The cutoff in the integral (12.81) from below at k — m, where m = 1/C is the reciprocal 
of the integral turbulence scale C, provides the IR regularization. Its precise form is 
unimportant; the sharp cutoff is the simplest choice for the practical calculations. 

The full problem is equivalent to the field theoretic model of the three fields 
$ = {if), if)\ v} with the action functional 

<?($) = ^t ( _ Vt + Xq9 2 _ XoTq)4) + ((^ _ ^2) + <j (v)) (2 .9) 

which is obtained from (12.31) by the replacement (12.61) and adding the term corresponding 
to the Gaussian averaging over the field v with the correlator (12. 8p : 

5(v) =~JdtJdxJ dx'vi^xjDf/^Vjifx'), (2.10) 



Ai(r) = D / -— - P fi (k) — f exp(ikr), k = |k| (2. 



where 



D^{r)^D^r-^ 



is the kernel of the inverse linear operation for the function -Dij(r) in (12.81) . 

In addition to (12. 5p . the Feynman diagrams for the model (I2.9P involve the 
propagator (w)o specified in (12 .7p . (I2.8P and the new vertex —ifj\vd)ifi. 

The role of the coupling constants in the ordinary perturbation theory is played by 
the two parameters 

u = gl ~ A 4 ~ d , = £>o/A ~ A 5 . (2.11) 

The last relations, following from the dimensionality considerations, (more precisely, see 
the next section) define the typical UV momentum scale A. By rescaling the fields, the 
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constant u>o can be placed in front of the interaction term —ip\vd)ip, which is more 
familiar for the field theory. We do not do it, however, in order not to spoil the natural 
form of the covariant derivative, and thus assign the factor wq to the propagator (vv)q. 

3. Canonical dimensions, UV divergences and the renormalization 

It is well known that the analysis of UV divergences is based on the analysis of canonical 
dimensions ("power counting"); see e.g. [HI [12]. Dynamic models of the type (12.91) . in 
contrast to static ones, have two independent scales: the time scale T and the length 
scale L. Thus the canonical dimension of some quantity F (a field or a parameter in the 
action functional) is completely characterized by two numbers, the frequency dimension 
d F and the momentum dimension dp, defined such that [F] ~ [T] - ^ [L] _d f\ These 
dimensions are found from the obvious normalization conditions 

d\ = —d* = 1, dk=d^ = 0, d k w = d k = 0, d^ = —d% = 1, 

and from the requirement that each term of the action functional be dimensionless (with 
respect to the momentum and frequency dimensions separately). Then, based on d k F and 
d F , one can introduce the total canonical dimension dp = dp + 2d F (in the free theory, 
d t oc d 2 ), which plays in the theory of renormalization of dynamical models the same 
role as the conventional (momentum) dimension does in static problems; see Chap. 5 of 

H2|. 

The dimensions for the model (12. 9p are given in table [H including renormalized 
parameters (without subscript "o"), which will be introduced later on. From tabled] 
or, equivalently, from the relations (12. lip it follows that the model is logarithmic (the 
both coupling constants go and il>o are simultaneously dimensionless) at d — 4 and 
£ = 0. Thus the UV divergences in the Green functions manifest themselves as poles in 
e = 4 — d, £ and, in general, their linear combinations. 

The total canonical dimension of an arbitrary 1-irreducible Green function T = 
i_i r is given by the relation [12] 

d r = 4 + 2c?r = d + 2 - N^, (3.1) 

where N$ = {N^, , N v } are the numbers of corresponding fields entering into the 
function T, and the summation over all types of the fields is implied. 
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The total dimension dr in logarithmic theory (that is, at e = £ = 0) is the formal 
index of the UV divergence 5r = dr\ e =^=o- Superficial UV divergences, whose removal 
requires counterterms, can be present only in those functions T for which Sr is a non- 
negative integer. From table [I] and (13. ip we find 

5 V = 6 - 2N* - 2N^ - N v . (3.2) 

We recall that in our model nonvanishing Green functions must involve the both fields 
ip and ip* (see the discussion in the preceding section), so in ( 13. 2ft it is sufficient to 
take > 1 and (simultaneously) iVt > 1. Straightforward analysis of the expression 
( 13. 2p then shows that superficial UV divergences can be present only in the following 
1-irreducible functions: 

(ip ip) (5 = 2) with the counterterms ip*d t ip, ip^d 2 ip, ip'ip, 
(ip^ipip) (5 = 0) with the counterterm ip ip 2 , 
(ip*ip*ip) (5 = 0) with the counterterm (ip*) 2 ip, 
(iP*iPv) (6 = 1), 

for which the counterterm necessarily reduces to the form ip*(vd)ip = —ip(vd)ip* owing 
to the transversality of the velocity field. All such terms are present in the action (12. 9p . 
so that our model appears multiplicatively renormalizable. 

The superficial divergence in the function (ip'ipvv) with 5 = and the counterterm 
ip'ipv 2 , allowed by the dimension, is in fact forbidden by the Galilean symmetry. 
Furthermore, the latter requires that the counterterms ifi'dtif) and ip^(vd)ip enter the 
renormalized action only in the form of the Lagrangian derivative if)* V^. 

Strictly speaking, the arguments based on the Galilean symmetry are applicable 
only to the velocity field governed by the Navier-Stokes equation, and generally become 
invalid for synthetic Gaussian velocity ensembles. It turns out, however, that for 
a Gaussian ensemble of the type (12.71) with vanishing correlation time the Galilean 
symmetry of the counterterms indeed takes place; see e.g. [27]. This issue, along with 
the consequences of the Galilean invariance for the renormalization, is discussed in the 
appendix A of Ref. [21] in detail. The proof given there is fully applicable to the model 
( 12. 9p . From the symmetry ( 12 Ah it also follows that the trilinear counterterms enter the 
renormalized action as the single combination (if) i ) 2 if) — ip*ip 2 . 

We thus conclude that the renormalized action can be written in the form 

S R (<$>) = ^ (-Z 1 V t + Z 2 \d 2 - Z 3 Xr) if) + 

+ Z^((iP*) 2 iP-iP*ip 2 )+S(v). (3.3) 

Here A, r and g are renormalized analogs of the bare parameters (with the subscripts "o") 
and fj, is the reference mass scale (additional arbitrary parameter of the renormalized 
theory). Since the last term <S(v) given by ( I2.10p is not renormalized, the amplitude Dq 
is expressed in renormalized parameters as 

D = w X = wXfi^. (3.4) 
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Expression (13. 3p can be obtained by the multiplicative renormalization of the fields 
ip — > ipZ^, ip^ — > ip'Z^i v — > vZ v and the parameters: 

A = XZ X , r = tZ t , g = gn e/2 Z g , w = w\i i Z w . (3.5) 

The renormalization constants in Eqs. (13. 3p and (13. 5p are related as follows: 

Z\ = Z,^Z^\ = Z v ZtyZ^\ Z 2 = ZjpZ^ Z\, Z% = Z^pZ^] Z\Z T , 

Zt± = Z^Z^ZgZx = Z^Z^ZgZx, 1 = Z W Z\. (3-6) 

Resolving these relations with respect to the renormalization constants of the fields and 
parameters gives 

Z v = 1, Z^p = Z^pj = Z-J , Z T = Z 3 Z 2 1 , 

Z\ = Z w x = Z 2 Z 1 1 , Zg = Z 4 Z 2 1 Z 1 ^ , (3-7) 

where the first equality is a consequence of the Galilean symmetry and the second - a 
consequence of the symmetry ( 12.41) . The first relation in the second line is a consequence 
of the absence of renormalization of the term <S(v) in (13.31) . For the coupling constant 
uo = #q introduced in ( 12. lip one has: 

u = u/j 6 Z u , Z u = Z 2 g . (3.8) 

The renormalization constants Z%—Z^ are calculated directly from the diagrams, then 
the constants in ( 13.51) are found from ( 13.71) and ( 13. 81) . 

The renormalization constants capture all the divergences at e,£ — > 0, so that the 
correlation functions of the renormalized model (13.31) have finite limits for e, £ = when 
expressed in renormalized parameters A, r and so on. In practical calculations, we used 
the minimal subtraction (MS) scheme, in which the renormalization constants have the 
forms Zi = 1+ only singularities in e and £, with the coefficients depending on the two 
completely dimensionless parameters — renormalized coupling constants u and w. To 
simplify the resulting expressions, it is convenient to pass to the new couplings, 

u — > u/16n 2 , w — > w/16n 2 ; (3.9) 

in what follows they will be denoted by the same symbols. 

Explicit calculation in the first order in u and w (one-loop approximation) gives 
rather simple results: 

Zi = l + T> Z * = 1 + l- 3 -ZF> ^3 = 1 + ^, Z 4 = l + H (3.10) 

4e 8e 4£ 2e e 

with the corrections of second order in u and w and higher. For w = one obtains (up 
to the notation) the well-known one-loop result for the model (j2.3P ; cf. e.g. [3]), while 
for u = the exact result for the rapid-change model is recovered; cf. e.g. [28J. 

4. RG functions and RG equations 

Let us recall an elementary derivation of the RG equations; detailed exposition can 
be found in monographs [TTJ [12] . The RG equations are written for the renormalized 
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Green functions G R = which differ from the original (unrenormalized) ones 

Q = ($■••$) only by normalization (due to rescaling of the fields) and choice of 
parameters, and therefore can equally be used for analyzing the critical behaviour. The 
relation Sn($,e,fj) = iS($,e ) between the functionals (12. 9p and (13. 3p results in the 
relations 

G(e ,...) = Z^Z*fG R (e } ^...). (4.1) 

between the Green functions. Here, as usual, and ALt are the numbers of 
corresponding fields entering into G (we recall that in our model Z v = 1); eo = 
{A , r , Mo, u) } is the full set of bare parameters and e = {A,r, u,w} are their 
renormalized counterparts; the dots stand for the other arguments (times/frequencies 
and coordinates/momenta). 

We use to denote the differential operation [id^ for fixed e and operate on both 
sides of the equation (14.11) with it. This gives the basic RG differential equation: 

{V RG + + } G R (e, n,...) = 0, (4.2) 

where T> R q is the operation T>^ expressed in the renormalized variables: 

V RG EE V„ + AA + (3 w d w - 7a Pa - lrV T . (4.3) 

Here we have written T> x = xd x for any variable x, and the anomalous dimensions 7 are 
defined as 

7_f ee In Zp for any quantity F, (4.4) 

and the (5 functions for the two dimensionless couplings u and w are 

(3 U ee V^u = u[-s - j u ], (3 W = V^w = w [-£ - j w ], (4.5) 

where the second equalities come from the definitions and the relations (13.51) . 

Equations (13.71) result in the following relations between the anomalous dimensions 
(S3]): 

7v = 7^t = 7i/2, 7a = -7u> = 72 - 7i, lv = 0, 

7r = 73 - 72, 7« = 2 74 - 272 - 7i- (4-6) 

The anomalous dimension corresponding to a given renormalization constant Zp is 
readily found from the relation 

If = (AA + AA) In Z F ~ - (eV u + £V W ) In Z F . (4.7) 

In the first relation, we used the definition (I4.4p . expression ( 14. 3 p for the operation 

in renormalized variables, and the fact that the Z's depend only on the two completely 

dimensionless coupling constants u and w. In the second (approximate) relation, we 

retained only the leading-order terms in the (3 functions ( 14. 51) . which is sufficient for the 

first-order approximation. The factors e and £ in ( 14.71) cancel the corresponding poles 

contained in the expressions (I3.10p for the constants Zp, which leads to the final UV 

finite expressions for the anomalous dimensions. This gives: 

u u 3w u 

7i = -4, 72 = - g + -j-, 7s = 74 = -u, (4.8) 
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and from (14. 6p one finally obtains: 



u u 3w 

7</> = 70t = - g , 7a = -J w = g + — , 

3u 3w 3u 3w , . 

7r = - T " T , 7« = " T - T , (4-9) 

with corrections of order w 2 , w 2 , uw and higher. 



5. Fixed points and scaling regimes 

It is well known that the long-time large-distance asymptotic behaviour of a 
renormalizable field theory is determined by the IR attractive fixed points of the 
corresponding RG equations. In general, coordinates of the possible fixed points are 
found from the requirement that all the (3 functions vanish. In the model (I2.9P the 
coordinates u#, w* are determined by the two equations 

Pu(u*, w*) = 0, /3 TO (u*, w*) = 0, (5.1) 

with the (3 functions given in (I4.5p . The type of a fixed point is determined by the 
matrix 

n = {tt ij = dfa/dg j }, (5.2) 

where fa is the full set of the j3 functions and gj = {g,w} is the full set of coupling 
constants. For an IR attractive fixed point the matrix Q is positive, that is, the real 
parts of all its eigenvalues are positive. 

From the definitions (14.51) and explicit expressions (14.91) for the anomalous 
dimensions we derive the following leading-order expressions for the /3 functions: 

fa = u{-e + 3u/2 + 3w/2), f3 w = w (-£ + u/8 + 3w/A). (5.3) 

From Eqs. (15.11) and ( 15.31) we can identify four different fixed points. For the three of 
them, the matrix Q appears triangular and its eigenvalues are given by the diagonal 
elements Q u = dfa/du and Q w = d[3 w /dw. 

1. Gaussian (free) fixed point: = u>* = 0; Q u = —e, Q w = — £ (all these expressions 
are exact). 

2. w* = (exact result to all orders), = 2e/3; Q u = e, Q w = — £ + e/12. 

In this regime, effects of the turbulent mixing are irrelevant in the leading-order IR 
asymptotic behaviour; the basic critical exponents are independent on £ and coincide 
to all orders with their counterparts for the "pure" DP class p]-[3]. However, the 
dependence on £ appears in the corrections to the leading-order behaviour, in particular, 
due to the correction exponent Q w . Although the expression for Q u is not exact (it has 
corrections of order e 2 and higher), the inequality Q u > is equivalent to e > within 
the e expansion. 
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3. = 0, = 4£/3 (exact); Q u = —e + 2£, Q w = £ (exact). 

In this regime, the nonlinearity (ip^) 2 ip — ip^ip 2 of the DP model is irrelevant, and 
we arrive at the rapid-change model of a passively advected scalar field ip [27] . For that 
model, the j3 function is given exactly by the one- loop approximation (see [28J ) , hence 
the exact results for w* and Q w . The dependence on e appears in the corrections, in 
particular, due to the correction exponent Q u . 

4. u* = 4(e — 2£)/5, w* = 2(12£ — e)/15. The eigenvalues of the matrix ( 15. 2ft have the 



20 V " v / v ' 

It is easily checked that they are both real for all e and £ (the expression in the square 
root is positive definite) and positive for e/12 < £ < e/2. 

This fixed point corresponds to a new nontrivial IR scaling regime (universality 
class), in which the nonlinearity of the DP model (12.31) and the turbulent mixing are 
simultaneously important; the corresponding critical exponents depend on the both RG 
expansion parameters e and £ and are calculated as double series in these parameters; 
see section [6l 

In figure [1] we show the regions of IR stability for all the fixed points in the e— £ 
plane, that is, the regions in which the eigenvalues of the matrix (15. 2p for the given fixed 
point are both positive. 

In the one-loop approximation (I5.3p . all the boundaries of the regions of stability are 
given by straight rays; there are neither gaps nor overlaps between the different regions. 
||| Such pattern is typical of the first-order approximations; cf. [7J [TOj [201 [211 I2H] for 
various models. The boundaries e < 0, £ < for point 1, e > for point 2 and £ > for 
point 3 are exact, while the other can be affected by the higher-order corrections: the 
boundaries will become curved and gaps or overlaps can appear between the different 
regions of IR stability. 

It is important that, for all these fixed points, the coordinates «*, w* are non- 
negative in the regions of their IR stability, in agreement with the physical meaning of 
the these parameters. It is also worth noting that the both boundaries for point 4 are 
determined by the same eigenvalue A~, which changes its sign at £ = e/2 and £ = e/12, 
while A + remains strictly positive in the entire region of stability. 

| For the first three fixed point this is obvious from the expressions for £l u ,w', for the point 4 this is 
quite unexpected at the first sight, but can be explained by the homogeneity of the expressions (|5.4|) 
in e and £, and is easily seen from the simple form ot the determinant, det f2 = (e — 2£)(12£ — e)/10. 



forms: 




(5.4) 
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£ = e/2 




£ = £r/12 



Figure 1. Regions of stability of the fixed points in the model 



6. Critical scaling and critical dimensions 

Existence of an IR attractive fixed point implies existence of scaling (self-similar) 
behaviour of the Green functions in the IR range. In this "critical scaling" all the 
"IR irrelevant" parameters (A, /i, g and w in our case) are fixed and the "IR relevant" 
parameters (coordinates/momenta, times/frequencies, r and the fields) are dilated. The 
critical dimensions Ap of the IR relevant quantities F are given by the relations 

A F = 4 + A^ + 7 *, A w = 2- 7 *, (6.1) 

with the normalization condition A^ = 1; see e.g. [12] for more detail. Here d k p W are 
the canonical dimensions of F, given in table [TJ and 7^ is the value of the anomalous 
dimension (14.41) at the fixed point: 7^ = 7^(11*, to*). This gives 

A^ = A^ t =rf/2 + 7 ;, A T = 2 + 7 ;. (6.2) 

Four fixed points of the model (I2.3P revealed in the preceding section correspond 
to four possible IR scaling regimes; for given e and £ only one of them is IR attractive 
and governs the IR behaviour. From the general expressions (16.11) . (16.21) and the explicit 
one-loop expressions (14.91) we find: 



1. Gaussian (free) fixed point; all the expressions are exact: 
= d/2, A T = A W = 2. 



(6.3) 
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2. Directed percolation (DP) regime; mixing irrelevant: 

A i , = 2- 7e/12, A r = 2 - e/4, A w = 2 - e/12. (6.4) 

These dimensions depend only on e, the corrections of order 0(e 2 ) and the references 
can be found in [3]. LSI The 0(e 3 ) calculation is in progress; see [30J. 

3. Obukhov-Kraichnan exactly soluble regime; all results exact (see e.g. [28] for the 
detailed discussion): 

A w = A r = 2 - e, = d/2. (6.5) 

4. New universality class (both mixing and DP interaction are relevant): 

A^ = 2 + (f - 3e)/5, A r = 2 - (e + 30/5, A w = 2 - f (exact). (6.6) 

The first two dimensions have nontrivial higher-order corrections in £ and £. The exact 
results for A^ in (16.51) and (16.61) follow from the general relations 7a = j w in (14. 6f) and 
A w = 2 — 7^ in (16.11) and the identity j w = £, which is a consequence of the fixed-point 
equation j3 w = with (3 W from (14.51) for any fixed point with ^ 0. 

Let us illustrate the consequences of these general scaling relations for the spreading 
of a cloud of the agent (or a cloud of "infected" particles) in the turbulent environment. 
The mean-square radius R(t) at time t > of a cloud of such particles, which started 
from the origin x' = at time t' = 0, is related with the linear response function in the 
time-coordinate representation as follows: 

**( t ) = /^ G ( t ,x), G(t ,x) = W ,^(0,0)), x = |x|, (6.7) 

see e.g. [3lJ. For the response function, the relations (16.3p - (l6.6p result in the following 
IR asymptotic expression: 

G(t,x)=x- 2 ^F(^-,^-), (6.8) 

with some scaling function F. Substituting (16. 8p into (16.71) gives (with the assumption 
that the integral converges) the desired scaling expression for the radius: 

where the scaling function / is related with F from (I6.8P as follows: 

f(z)= J dxx 2 - 2A *F(x,4 

Directly at the critical point (assuming that the function / is finite at r = 0) one obtains 
from (16. 9p the power law for the radius: 

R 2 (t) OC ^+2-2A„)/A u = t (2-2 7 ;)/(2- 7 *). (OQ) 



S The conventional critical exponents used e.g. in [3] are related to the critical dimensions from (|6.4 
as z = A w , 1/v = A r , d + rj — 2 . 
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in the second equality we used the relations (16.1 j) and (16. 2p . For the Gaussian fixed 
point (16. 3p the usual "1/2 law" R(t) oc t 1 / 2 for the ordinary random walk is recovered. 
For the fixed point (16.51) . where the DP nonlinearity is irrelevant, one obtains the exact 
result R(t) oc t l ^ 2 ~^ . For the most realistic (Kolmogorov) value £ = 4/3 this gives 
R(t) oc t 3//2 in agreement with Richardson's "4/3 law" dR 2 /dt oc i? 4//3 for a passively 
advected scalar impurity; see [31]. For the other two fixed points the exponents in (16. 9p . 
(I6.10P are given by infinite series in e (point 2) or e and £ (point 4); the first-order 
approximations are easily obtained from (16.41) and (16. 6p . 

7. Conclusion 

We studied effects of turbulent stirring and mixing on the reaction-diffusion system in 
which the spreading of an agent (e.g., infectious disease) occurs, as an example of critical 
behaviour in an nonequilibrium system near its transition between the absorbing and 
fluctuating phases. We coupled two paradigmatic models: the so-called simple epidemic 
process, also known as Gribov's process and equivalent to Reggeon field theory (I2.3p . 
directed bond percolation process, Schlogl's first reaction and so on p] [3], and the 
Obukhov-Kraichnan's rapid-change model (I2.8P for the advecting turbulent velocity 
field [27]. The full problem can be reformulated as a multiplicatively renormalizable 
field theoretic model (12.91) . which allows one to apply the field theoretic renormalization 
group [TU [12] to the analysis of its critical behaviour. 

We showed that, depending on the relation between the spatial dimension d and the 
exponent £ that comes from the Obukhov-Kraichnan's ensemble, the model exhibits four 
different critical regimes, associated with four fixed points of the RG equations. Three 
fixed points correspond to known regimes: Gaussian fixed point (ordinary diffusion or 
random walk), directed percolation process with no advection (DP class), and passively 
advected scalar field (reaction processes, described by Gribov's nonlinearity, appear 
irrelevant). The fourth point reveals existence of a new nonequilibrium universality class, 
in which both the reaction and the turbulent mixing are relevant; the corresponding 
critical exponents are calculated to the leading order (one-loop approximation) of the 
double expansion in £ and e = 4 — d. 

Judging naively from the dimensions of the coupling constants (12. lip one could 
expect that the latter regime must take place when £ and e are simultaneously positive, 
but the careful RG analysis shows that the region of IR stability of the corresponding 
fixed point is much narrower (in the one-loop level it shrinks to the sector e/12 < £ < 
e/2). Of course, this is not a big surprise: similar phenomenon was established a long 
ago for the interplay between the long-range and short-range spin-spin interactions in 
models of equilibrium critical behaviour [32] . Later it was also encountered in the long- 
range versions of various models with nonequilibrium transitions [7J [33] and in models 
of a passive scalar advected by the velocity field with finite correlation time [29] . 

In our case, this effect leads to interesting prediction: in contrary to what could 
be naively anticipated, the most realistic spatial dimensions d = 2 or 3 and the 
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Kolmogorov's exponent £ = 4/3 for the fully developed turbulence correspond not to the 
most nontrivial new regime, but to the passive-scalar fixed point: the nonlinearity of the 
Reggeon model is irrelevant and the spreading of the agent is completely determined by 
the turbulent transfer. In particular, the time spreading or a cloud of infected particles 
(or of the agent) behaves accordingly to the power law R(t) oc t 1 ^ 2- ^, which is the 
proper generalization of Richardson's law to the case of arbitrary exponent £ in the 
velocity correlator (12. 8p . 

Further investigation should take into account anisotropy of the experimental set- 
up, compressibility, non-Gaussian character and finite correlation time of the advecting 
velocity field, effects of memory (immunization); interaction of the order parameter with 
other relevant degrees of freedom (mode- mode coupling), feedback of the reaction on 
the dynamics of the velocity and so on. This work is already in progress. 

Acknowledgments 

The authors thank Loran Adzhemyan, Michal Hnatich, Juha Honkonen and Mikhail 
Nalimov for helpful discussions. One of the authors (NVA) thanks H-K Janssen for 
bringing his attention to Refs. [TOj [33] . The work was supported in part by the Russian 
Foundation for Fundamental Research (grant No 08-02-00 125a), the Russian National 
Program (grant No 2.1.1.1112) and the program "Russian Scientific Schools" (grant 
No 5538.2006.2). 

Appendix A. Calculation of the renormalization constants 

In this section we derive the first-order results (13.101) for the renormalization constants. 
Of course, the one-loop calculation is rather simple and can be accomplished in a few 
different ways, but we will discuss it for completeness and in order to mention some 
interesting subtleties specific of the model (12.91) . 

The renormalization constants can be found from the requirement that the Green 
functions of the renormalized model fl 3 . 3 j) . when expressed in renormalized variables, be 
UV finite (in our case, be finite at e — > 0, £ —>■ 0). Owing to the symmetry i \2A\\ and to 
the Galilean invariance, the full set of constants Z\-Z^ in our model can be found from 
just two 1-irreducible functions: (-^Vh-ir aud (V^V* V) 1 - 1 * (or, if desired, {ip^ipip)^^). 
In the renormalized model, the corresponding one-loop approximations have the forms 

(^V)i-ir = i^Zi + Xp 2 Z 2 + XtZ 3 - i 

(A.l) 

(which is the Dyson equation for the exact function (ip^ip), hence the minus signs 
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in front of the diagrams) and 



(VVV)l-ir = ^4 + 2 / \+2/\+^/\ (A.2 




The first two diagrams in (1A.2D have in fact two different forms, related by the mirror 
reflection, but they give equal contributions to the renormalization constants and are 
accounted by the factors of 2. Owing to the Galilean symmetry, the constant Z\ can 
also be found from the function 




(^Vv)i-ir = -ip^i + / \ + / \ • (A.3) 



The solid lines in (jA.lj) — (jA.3j) denote the propagator (V>V^)o from (12. 5ft (the arrow is 



directed from ip to i[)>) and the wavy lines correspond to (vv)o from (12. 8p . All the 
diagrammatic elements should be expressed in renormalized variables using the relations 
(13.3p - (13.6D . In the one- loop approximation, the Z's in the bare terms of I1A.1I) . (1A.2I) 
should be taken in the first order in u = g 2 and w, while in the diagrams they should 
simply be replaced with unities, Zi — > 1. Thus the passage to renormalized variables in 
the diagrams is achieved by the simple substitutions A — > A, r — ■> r, g — > g/// 2 and 
Wo — > W[£. 

The second and third diagrams in ( 1A.2I) and the second diagram in ( 1A.3I) appear 
UV finite and therefore give no contribution to the renormalization constants. Indeed, 
due to the transversality of the velocity field, the derivative in the vertex —ip^(vd)ip 
can also be moved onto the field if}' using integration by parts: —ip^ (vd)if> = if)(vd)if>\ 
Thus in any diagram involving n external vertices of this type, the factor p n with n 
external momenta p will be taken outside the corresponding integrals. This reduces the 
dimension of the integrand by n units and can make it UV convergent. In the case at 
hand, this proves the UV finiteness of the three diagrams mentioned above: for all of 
them n = 2 while the formal index of divergence is or 1. 

What is more, since the propagator (12.51) is retarded and (12.81) contains the 5 
function in time, the second diagrams both in flA.21) and (1A.3|) contain self-contracted 



"circuits" of the step functions in time and therefore vanish identically. This argument, 
however, does not apply to the second diagram in (1A.1I) . which requires a more careful 
treatment; cf. [2lT [28] . The analytic expression for that diagram has the form 

f du f dk DoPijjk) 1 
PiPj J (2tt) J k>m (2nY k d+ Z -ico + a(p - k) ' ^ 
where the prefactor comes from the vertices, the IR cutoff and the first cofactor in the 
integrand come from the propagator (12.81) and the second cofactor with cr(k) = A(/c 2 + r) 
is the propagator (12. 5p . The expression (1A.4|) is independent of the external frequency. 
Integration over uo involves the indeterminacy 

7^—^ xA= d ^ (A- 5 ) 

(27r) — \U) + er(p — k) 
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where 6(0) is the step function at the origin. This reflects the details of the velocity 
statistics lost in the white-noise limit; see the discussion in [27] . In the case at hand, 
the S function in ( 12. 7p should be understood as the limit of a narrow function which 
is necessarily symmetric in t, t', because one deals with a pair correlator. Thus the 
indeterminacy in (1A.5I) must be unambiguously resolved as half the sum of the limits: 
8(0) = 1/2. Then the remaining integral over k in (IA.4I) appears independent of p and 
r and is easily calculated: 

] k>m (2*y d J k>m (2K)*k*< ijm (2ttY d$ ' 

where Sd = 2n d / 2 /T(d/2) with Euler's T function is the surface area of the unit sphere 
in the <i-dimensional space. Collecting all the factors and setting D = XoWo = Xw^ 
gives the final result for the diagram: 

- AV(A0») € = -Vip^ + UV finite part. (A.6) 

In the last equality in (1A.6j) . the UV divergent part of the diagram is selected (the 
replacements (n/mp — > 1 and d = 4 — e — > 4 are made); it contains a first-order pole in 
£. The expression ( 1A.6I) as a whole is proportional to p 2 , so that it gives a contribution 
only to the renormalization constant Z 2 in (lA.lj) . 

The remaining two diagrams in ( lA.lj) . ( 1A.2j) do not involve the velocity correlator; 
they are independent of £ and contain only poles in e. Although the calculation of these 
diagrams is discussed in [3] within the context of the model (12.31) . we will sketch an 
alternative calculation here, mainly in order to present a reference formula, which can 
be interesting in itself. 

The key point is as follows: the convolution of two functions of the form 

F(a,a,r) = (-icu a + k 2 + r)' a , r>0 (A.7) 
is a function of the same form: 

F(ai, a u Ti) * F(a 2 , a 2 , r 2 ) = K(a 1 , a 2 ; a t , a 2 ) F(a 3 , a 3 , r 3 ) (A. 8) 
if ai and a 2 have the same sign and zero otherwise. Here 

a 3 = ai + a 2 , a 3 — at + a 2 — d/2 — 1, r 3 = a 3 {r\/ai + r 2 /a 2 ); (A. 9) 
the coefficient is independent on the T\ )2 and has the form: 

jw \ _ d/2-ai d/2-a 2 ai +a 2 -d-l F(a 3 ) 

A («!, a 2 , oi, a 2 ) - a l a 2 a 3 ^- / -^^. 

For the "massless" case ri j2 = this formula was proposed and used in the three- 
loop calculation of the critical exponent z in the so-called model A of critical dynamics 
in [33]; see also [2U [35] for other applications. For the general case (ti j2 > 0) equation 
(IA.8j) can be obtained from the following observation: In the time-space representation 
the function ( 1A.7I) takes on the form 

^/ \ 0(t sign(a)) a d / 2 ~ a „ H/9 -, f ax 2 .1 . . 

n«,*,r) - 1 exp i it - Tt/a s ' (A - 10) 
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and the product of such functions (which corresponds to the convolution of their Fourier 
transforms) is obviously a function of the same form. Note that, when a\ and a<i have 
different signs, the convolution in ( 1A.8I) corresponds to the product of a retarded and 
an advanced functions of the form (lA.lOj) and therefore vanishes. 

Let us turn to the remaining diagrams in (1A.1I) . ( 1A.2I) . With no loss of generality, 
one can set A = 1 (the dependence on A can be restored in the final answers 
using the dimensionality considerations). Then the frequency-momentum integral that 
corresponds to the first diagram in ( 1A.1I) . is equal to the convolution ( 1A.8I) . in which 
a lj2 = 1, ai j2 = 1, Ti )2 = t. From (1A.9I) one obtains a 3 = 2, a 3 = — 1 + e/2, r 3 = At, so 
that the left-hand side of (1A.8I) is proportional to 

r(-l + e/2)(-2iw + k 2 + At) 1 -" 12 = ~(-2iu + k 2 + At) + UV finite part 

(A.ll) 

(here and below we omit uninteresting finite numerical factors, powers of 7r and so on). 

The first diagram in ( 1A.2I) is logarithmically divergent. According to the general 
statements of the renormalization theory, its divergent part in the MS scheme does not 
depend on the specific choice of the external momenta, frequencies, and "masses" like 
r, provided this choice guarantees IR convergence of the corresponding integral; see e.g. 
[TT| [12] . Thus we can set the external frequency and momentum flowing into the right 
lower vertex equal to zero. (We could also set r = 0, but this is not necessary). Then 
the integral that to the diagram in the frequency-momentum representation, becomes 
equal to the convolution (1A.8I) with = 1, ci\ = 1, «2 = 2, T\,i = r. From ( 1A.9I) it 
follows a 3 = 2, a 3 = e/2, r 3 = At, and the left-hand side of (1A.8|) is proportional to 

r(e/2)(-2iw + k 2 + At)- £/2 = ~ + UV finite part. (A.12) 

As expected, the expressions ( 1A.11I) and (IA.12I) contain first-order poles in e. Their 
pole parts are polynomials in frequencies, momenta, and "masses," so that they can 
be cancelled in expressions ( lA.lj) . (1A.2I) by the proper choice of the renormalization 
constants Z\—Z±. Taking into account all the factors (signs, symmetry coefficients, 
factors ±g from the vertices) and the replacement (13. 9ft gives the results announced in 
( I3.10p . It remains to note that the calculation of the 1-irreducible function (1A.3|) indeed 
gives the same result ( 13.101) for the constant Z±, in agreement with general consequences 
of the Galilean invariance of our model. 
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